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The ultimate goal of the observation of nonthermal emission from astrophysical sources is to un¬ 
derstand the underlying particle acceleration and evolution processes, and few tools are publicly 
available to infer the particle distribution properties from the observed photon spectra from X- 
ray to VHE gamma rays. Here I present naima, an open source Python package that provides 
models for nonthermal radiative emission from homogeneous distribution of relativistic electrons 
and protons. Contributions from synchrotron, inverse Compton, nonthermal bremsstrahlung, and 
neutral-pion decay can be computed for a series of functional shapes of the particle energy dis¬ 
tributions, with the possibility of using user-defined particle distribution functions. In addition, 
naima provides a set of functions that allow to use these models to fit observed nonthermal spectra 
through an MCMC procedure, obtaining probability distribution functions for the particle distri¬ 
bution parameters. Here I present the models and methods available in naima and an example of 
their application to the understanding of a galactic nonthermal source. 

naima’s documentation, including how to install the package, is available at 
http://naima.readthedocs.org. 
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1. Introduction 

Over the past few years, there has been a wealth of facilities that allow for an unprecedented 
level of sensitivity in the high-energy and very-high-energy gamma-ray bands, most notably Fermi- 
LAT and the Cherenkov telescope arrays H.E.S.S., MAGIC, and VERITAS. These energy bands 
probe the highest energies at which particles are known to be accelerated, and provide valuable 
insight into the acceleration and energy loss processes of a growing population of galactic and 
extragalactic sources. 

A first step in understanding the mechanisms through which a given source has accelerated the 
particles that are responsible for the emission we detect is to characterize the present age particle 
energy distribution. Eor many sources, such as pulsar wind nebulae emitting gamma-rays through 
inverse Compton scattering on the Cosmic Microwave Background (CMB), the derivation of the 
particle distribution can be done in a model-independent way. To date, there are few public codes 
that allow this sort of analysis, and in many papers the discussion of observed spectra is carried 
out based on modelling code that is not published along with the paper, therefore hampering the 
reproducibility of the research. 

Here I present naima, an open source Python package developed with the aim of providing 
proven methods for computing the nonthermal radiative output of relativistic particle distributions. 
Being open source makes it easy for anyone to check the algorithms used in the calculations, and 
of improving and extending the package. Python is increasingly gaining ground as the language 
of choice for astronomical research, and makes it easy to combine existing libraries such as naima 
into a given analysis. 

There are two main components of the package: a set of nonthermal radiative models, and a 
set of utility functions that make it easier to fit a given model to observed spectral data. 

2. Features 

The radiative models are implemented in a modular way that allows to select a functional shape 
for the particle distribution. There are several functions included for this purpose, and the user can 
as well define their own following the instructions in the documentation, therefore allowing to 
implement any type of particle cooling, escape, or acceleration physics before computing its radia¬ 
tive output with naima. The radiative models currently available in naimaare synchrotron, inverse 
Compton, nonthermal bremsstrahlung, and neutral pion decay. See Section 3 for details on the im¬ 
plementation. In addition, naima includes a set of wrappers around these models that allow them 
to be used within the sherpa spectral analysis package [1] in the module naima. sherpa_models. 

A set of fitting utilities are provided in naima with the goal of infering the properties of the 
parent particle distribution that gives rise to an observed nonthermal spectrum. These use Markov 
Chain Monte Carlo (MCMC) sampling to obtain both the maximum likelihood parameters as well 
as their uncertainties in a single run. Details on the implementation can be found in Section 4. 
Several plotting functions are available as well to analyse the results of the MCMC run. 

naima uses the astropy package [2] extensively, most notably through the use of its physical 
unit module astropy. units. This allows users to define the input spectra and parameters in their 
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preferred units, and naima will be able to convert them as needed in its internal calculations, with 
the added benefit of ensuring that the algorithms in naimaare dimensionally correct. 

3. Implementation of radiation models 

3.1 Synchrotron 

Synchrotron radiation is produced by all charged particles in the presence of magnetic fields, 
and is ubiquitous in the emitted spectrum of leptonic sources. A full description and derivation of 
its properties can be found in [3]. 

The naima.models.Synchrotron class implements the parametrization of the emissivity 
function of synchrotron radiation in random magnetic fields presented by [4, appendix D]. This 
parametrization is particularly useful as it avoids using special functions, and achieves an accuracy 
of 0.2 entire range of emission energy. 

3.2 Inverse Compton 

The inverse Compton (IC) scattering of soft photons by relativistic electrons is the main 
gamma-ray production channel for electron populations [3]. Often, the seed photon field will 
be a blackbody or a diluted blackbody, and the calculation of IC must be done taking this into 
account, naima implements the analytical approximations to IC upscattering of blackbody radia¬ 
tion developed by [5]. These have the advantage of being computationally cheap compared to a 
numerical integration over the spectrum of the blackbody, and remain accurate within one percent 
over a wide range of energies. Both the isotropic IC and anisotropic IC approximations are avail¬ 
able in naima. Note that even environments with complex broadband background radiation can be 
effectively modelled by a combination of diluted blackbodies. 

The implementation in naima allows to specify which blackbody seed photon fields to use in 
the calculation, and provides the three dominant galactic photon fields at the location of the Solar 
System through the CMB (Cosmic Microwave Background), FIR (far-infrared dust emission), and 
NIR (near-infrared stellar emission) keywords. 

3.3 Nonthermal Bremsstrahlung 

Nonthermal bremsstrahlung radiation arises when a population of relativistic particles interact 
with a thermal particle population. For the computation of the bremsstrahlung emission spectrum, 
the Bremsstrahlung class implements the approximation of [6] to the exact cross-section pre¬ 
sented by [7]. Electron-electron bremsstrahlung is implemented for the complete energy range, 
whereas electron-ion bremsstrahlung is at the moment only available for photon energies above 
10 MeV. The normalization of the emission, and importance of the electron-electron versus the 
electron-ion channels, can be selected in the class, with the default values assuming a fully ionised 
target medium with solar abundances. 

3.4 Pion Decay 

The main gamma-ray production for relativistic protons are p-p interactions followed by pion 
decay, which results in a photon with Ey > 100 MeV. Until recently, the only parametrizations 
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available for the integral cross-section and photon emission spectra were either only applicable 
to limited energy ranges [8, > 0.1 TeV], or were given as extensive numerical tables [9]. By 

considering Monte Carlo results and a compilation of accelerator data on p-p interactions, [10] 
were able to extend the parametrization of [8] down to the energy threshold for pion production. 
The PionDecay class uses an implementation of the formulae presented in their paper, and gives 
the choice of which high-energy model to use (from the parametrization to the different Monte 
Carlo results) through the hiEmodel parameter. 


4. MCMC sampling 


The following will briefly describe the implementation of spectral fitting in naima, and a full 
explanation of MCMC and the sampling algorithm can be found in [11], and in the documentation 
of emcee, the package used for MCMC sampling [12]. 

The measurements and uncertainties in the observed spectrum are assumed to be correct, Gaus¬ 
sian, and independent (note that this is unlikely to be the case, see Section 6 on how this might be 
tackled in the future). Under this assumption, the likelihood of observed data given the spectral 
model S (p; E), for a parameter vector p, is 


^=n 
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where {Fi,cri) are the flux measurement and uncertainty at an energy Ei over N spectral mea¬ 
surements. Taking the logarithm. 
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Given that the MCMC procedure will sample the areas of the distribution with maximum 
value of the objective function, it is useful to define the objective function as the log-likelihood 
disregarding constant factors: 
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The InX function in this assumption can be related to the;^2 parameter as;^2 - -21nX, so that 
maximization of the log-likelihood is equivalent to a minimization of;^2 

In addition to the likelihood from the observed spectral points, a prior likelihood factor should 
be considered for all parameters. This prior likelihood encodes our prior knowledge of the prob¬ 
ability distribution of a given model parameter. The combination of the prior and data likelihood 
functions is passed onto the emcee sampler function, and the MCMC run is started, emcee uses an 
affine-invariant MCMC sampler that has the advantage of being able to sample complex parameter 
spaces without any tuning required. In addition, having multiple simultaneous walkers improves 
the efficiency of the sampling and reduces the number of computationally-expensive likelihood 
calls required. 
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Figure 1: Diagnostic plot produced by naima.plot_chain for the sampling of the cutoff energy in the 
particle distribution. The top-left panel shows the traces for the 32 walkers in gray, with three of them 
highlighted in red, that can be used to estimate whether the sampling has stabilized around the maximum 
likelihood parameters. The right panel shows the posterior distribution of the parameter, along with the 
median (dashed line) and an estimate of the Icr confidence interval (gray band). On the bottom-left, statistics 
of the parameter distribution, including a median with uncertainties based on the 16th and 84th percentiles, 
is shown. Note that the parameter is sampled in logarithmic space in this example, and given the label 
log 10(cutoff) naima can identify this and convert it to its value and uncertainty in TeV, resulting in the 
estimate £’cutoff = 11^-20 


5. Example analysis: hadronic emission from RXJ1713-3946 

As an example of nonthermal spectral analysis with naima, here we will show the results 
of inferring the particle distribution parameters of a hadronic population from a spectrum of the 
shell-type supernova remnant RXJ1713-3946 obtained with the H.E.S.S. Cherenkov telescope 
array [13]. The photon spectrum in the 0.3-100TeV energy range is well characterized by a 
power-law with an exponential cutoff, so we will use a similar function for the particle distri¬ 
bution. Even though the leptonic or hadronic origin of the gamma-ray enfission from RX J1713- 
3946 is still under debate [14, and references therein], for the sake of example here we will as¬ 
sume an neutral pion decay origin of the VHE gamma-ray enfission. Several other examples, 
with the full source code needed to reproduce them, are available in the naimadocumentation 
http://naima.readthedocs.org. 

The first step is the definition of the radiative model to be fit, and in this case we will use 
the ExponentialCutof fPowerLaw particle distribution function and PionDecay radiative model. 
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Figure 2: Distribution of the free model parameters, norm is the particle distribution normalization at 5 TeV 
in units of TeV“\ index is the power law index, and log 10 (cutoff) is the decimal logarithm of the cutoff 
energy in units of TeV during the MCMC sampling run. The parameter vector with the highest likelihood 
found during the run is indicated by the red cross. 

The spectrum from [13] can be read as a table with astropy. io. ascii, and passed onto the naima 
MCMC sampling functions along with the defined model function. We run the sampling with 32 
simultaneous walkers (or sampling chains), for 100 steps of burn-in, and 400 steps of run that are 
saved for later analysis. The resulting sampled chains can be analysed with a set of naima functions 
that plot diagnostic and results plots. Such a chain diagnostic produced with naima. plot_chain 
can be seen for the energy of the exponential cutoff in the particle distribution model in Figure 1 . 
Figure 2 shows a corner plot, which plots the distribution for all parameters against each other. It 
can be seen that there is a large positive covariance between the distributions of the particle index 
and cutoff parameters. Finally, Figure 3 shows a comparison between the observed spectrum and 
the computed model, including the spread of the spectrum derived from 100 random parameter 
vectors taken from the MCMC chain, and an inset with the inferred particle energy distribution. 

The result of the MCMC run provides estimates of the parameters of the particle distribution 
function. For this spectrum and radiative model, the particle index is constrained to ^ = 1.90 ± 
0.08, and the cutoff energy is MCMC parameter chain can also be used 

to compute the distribution of quantities derived from the particle spectrum, such as total energy 
content in protons, which for RXJ1713-3946 is constrained to Wp(Ep > ITeV) = (5.74;^^ 2 ^) x 
erg, where nu is the target density and <iikpc is the distance to the source in kpc. 

6. Limitations and future development 

The main limitation of the approach used by naima for spectral fitting is the assumption of 
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Figure 3: H.E.S.S. spectrum of RXJ1713-3946 [13], computed spectrum from a hadronic model, and 
residuals of the maximum likelihood model (bottom panel). The thick black line indicates the maximum 
likelihood spectrum, and the gray lines are 100 samplings of the posterior distribution of the model parameter 
vector. The inset shows the energy distribution of the proton population in erg versus the proton energy in 
TeV. 


uncorrelated, gaussian errors. Even though this may be incorrect for many spectra, mostly when 
considering fine structure, it is often the only approach possible when simultaneously fitting pub¬ 
lished spectra from radio to VHE gamma-rays. When instrument response functions are available, 
a way to avoid this assumption is to use the sherpa models in naima. sherpa_models. 

Future development of the package will focus on the addition of simple particle cooling func¬ 
tions (more complex physics, such as time-dependent particle evolution, should be done on a 
case-by-case basis), and use of naima radiative models in gammapy [15], a Python package for 
gamma-ray data analysis. 
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